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1 Introduction 

The engineering uses of rubber have expanded well beyond traditional products such as tires and seals. 
Today rubber components can be found in a diverse set of constructs including engine mounts, building 
foundations, belts, and fenders (see [12], [21]). Increasingly, the applications of rubber are becoming more 
sophisticated, as exemplified by the use of rubber bearings in bridges which allow for thermal expansions of 
the deck without placing excessive loads on the bridge supports (see [14]). 

In current engineering applications, rubber or elastomer composites are typically filled with inactive 
particles such as carbon black or silica. If active fillers were used, such as piezoelectric, magnetic, or 
conductive particles, the resulting controllable elastomer could be used in products such as active vibration 
suppression devices (e.g., see [16], [17]). As these new materials are developed, the role of design will increase 
in both complexity and importance. In particular, the capability to predict the dynamic mechanical response 
of the components will become increasingly valuable. 

The many desirable characteristics of rubber as a design component, which include the ability to undergo 
large elastic deformations and provide significant damping with near incompressibility, are also contributing 
factors to the complications arising in the process of formulating models. Damping is highly complex, 
and the strain history, rate of loading, environmental temperature, and amount and type of filler affect 
the mechanical response in a nontrivial manner. Additionally, many elastomers exhibit strong hysteresis 
characteristics similar to those found in shape memory alloys and piezoceramic actuators. 

Lord Corporation, a company based in Cary, NC, develops and manufactures a wide variety of damping 
devices, including many types of elastomer-based dampers. In 1994, applied mathematicians with the Center 
for Research in Scientific Computation (CRSC) at North Carolina State University began collaborations with 
scientists at Lord Corporation to develop a high-resolution model of the dynamical behavior of elastomers. 
The ultimate goal of this project is to construct a dynamic model and accompanying design tools that will 
assist Lord in the development of higher quality damping devices. The original team members included H.T. 
Banks, N.J. Lybeck, Y. Zhang and N. Medhin from the CRSC, and M.J. Gaitens, B.C. Munoz, and L.C. 
Yanyo from Lord Corporation. 

Our team has concentrated on two types of deformations: simple extension and simple shear. Since 
most complex motion is believed to be a combination of these two types of motion (see [9], [22]), a basic 
understanding of each is an important precursor to a full understanding of complex deformations. Here we 
will give a summary of our efforts to date on simple extension. 
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The first two years of the team’s collaboration included an extensive study of existing elastomer models. In 
spite of the many complex issues relating to elastomer dynamics, researchers have made substantial progress 
in developing tools for such models (see [13],[26],[27] for basic texts). These models are predominantly 
phenomenological, based on strain energy function (SEF) and finite strain (FS) theories. SEF theories contain 
information about the elastic properties of elastomers, but do not describe either damping or hysteresis, and 
hence are typically used for static finite element analysis (see [9]). The CRSC/Lord team worked, both 
theoretically and computationally, with several of these models. Computations based on the existing models 
led to the development of more general nonlinear models. Dynamic and quasi-static experiments were 
designed and carried out to test the performance of the models on different types of elastomers. Details on 
these models and experiments will be outlined in Section 2. 

In working with the existing and generalized nonlinear models, the team recognized the need for a model 
that could incorporate the effects of hysteresis and damping. Additional experiments were conducted, and 
the quasi-static case was considered for developing a nonlinear, hysteretic stress-strain constitutive law. G.A. 
Pinter and L.K. Potter of the CRSC joined the team in 1996 as hysteresis was being incorporated into both 
the quasi-static and dynamic models. The formulation of the dynamic model led to the question of well- 
posedness, which was addressed and resolved. Finally, additional computations were performed to validate 
the accuracy of the model on several types of elastomers. These results will be presented in Section 3. 

2 Nonlinear extension models, experiments and results 

2.1 Neo-Hookean extension models 

The strain energy function (SEF) and finite strain (FS) models are often used to develop a dynamic model 
for the behavior of elastomers. These models are designed for materials that are usually assumed to be 
isotropic and incompressible; the simplest of these are known as neo-Hookean materials. 

The SEF material models use the principal extension ratios A i to represent the deformed length of 
unit vectors parallel to the principal axes (the axes of zero shear stress). The SEF models of Mooney 
and Rivlin are based on Rivlin’s proposal [23] that the SEF should depend only on the strain invariants 
II = A^ + X'i + A§, In = AoAg + A^A| + A^Ao and I3 — A^AoA^. For example, the Mooney SEF is given 
by U = ( '• (/• — 3) + Cn[In — 3), or more generally, the modified expression U = ( '•(/• — 3) + f(In — 3), 
where / has certain qualitative properties. This expression is most appropriate for components where the 
rubber is not tightly confined and where the assumption of absolute incompressibility (implying A 1 A 2 A 3 = 1 
or I 3 = 1) is a reasonable approximation. The more general Rivlin SEF U = ^ — 3)*(/2 — 3) J 

and its generalization for near incompressibility (see [9]) allow higher order dependence of the SEF on the 
invariants. The works of Ogden [22], as well as Valanis and Landel [9], represent an important departure 
from Rivlin’s proposal, in that their formulations are based on SEFs that depend only on the extension 
ratios. 

The finite strain elastic theory of Rivlin [27],[23] is developed with a generalized Hooke’s law in an analogy 
to infinitesimal strain elasticity, but requires no “small deformation” assumption and includes higher order 
exact terms in its formulation. Moreover, finite stresses are defined relative to the deformed body and hence 
are the “true stresses”, as opposed to the “nominal” or “engineering” stresses (relative to the undeformed 
body) one usually encounters in the infinitesimal linear elasticity used with metals. This Eulerian measure of 
strain is an important feature of any development of models for use in analytical/comp utaiional/expen mental 
investigations of rubber-like material bodies. The finite strain elasticity of Rivlin can be directly related to 
the strain energy function formulations through equations relating the finite strains e xx , e yy , e, z to the 
extension ratios A 1; A 2 , A 3 used to define the SEF. 

The finite strain approach can be put in a somewhat more general perspective in the context of classical 
modeling of elastic solids and fluids (filled elastomers do not fit exactly into either category, although in 
the absence of cross links, elastomers are highly viscous fluids). In classical approaches one frequently 
encounters an Eulerian formulation in dynamics of fluids where large deformations or displacements are 
common, whereas a Lagrangian formulation is employed for solids undergoing small elastic deformations. 


2 



In both formulations, momentum balance laws along with constitutive laws relating stress and strain are 
employed to develop theories of dynamics (see [22],[20]). In the general Lagrangian formulations, quantities 
(such as stress, strain) are defined relative to an original or reference configuration Bo of the body or structure 
in terms of a fixed coordinate system X = ]Ji ,}. For an Eulerian formulation one defines quantities relative 
to a “current” or deformed configuration B with coordinates x — {aq} relative to the deformed configuration. 

A fundamental role in discussing the relationship between these formulations is a “configuration” map 
x — x(A') or “motion” x(t) = x(t, X) if the deformations are changing in time from an original configuration 
x(0, X) = X. The deformations (in the usual elasticity terminology) are then given by 

u(t, X) = x(t, X) - x(0, X) = x(t, X) - X. 


by 


The configuration gradient (also called the “deformation gradient” in an unfortunate misnomer) is defined 

A _ dx{X) _ dx_ 
dX dX 


and is used to define the right (T t T) and left (AA T ) Cauchy-Green “deformation tensors” and the usual 
(in elasticity theory) Green Si. Venant strain 


£ = i(T T T - I) = ^(D t D + D + D t ), 


where D = J^|- = A — I. These definitions, along with momentum balance laws, can be used to write 
dynamic equations in either the Lagrangian (reference) or Eulerian (current) coordinate systems. If one uses 
the Eulerian system, stresses are given in terms of the Cauchy or true stress T and constitutive (stress- 
strain) laws are expressed in terms of current coordinates x. For computational and experimental purposes, 
however, it is often more desirable to use a Lagrangian coordinate system and then the Cauchy or true stress 
must be converted to an expression for the nominal or engineering stress S = JA~ 1 T , where J = det. A and 
S T is called the First Piola-Kirchoff stress tensor. 

The SEF or finite strain formulations can be used along with standard material independent force and 
moment balance derivations (the Timoshenko theory [25],[11]) as the basis of a dynamic model for extension. 
Here we take a simple example: an isotropic, incompressible (A 1 A 2 A 3 = 1) elastomeric rod with a tip mass 
under simple elongation with a finite applied stress in the principal axis direction x\ = x (see Figure 1). 
Following standard convention, we use lower case letters to denote the Lagrangian coordinates. The finite 
stress theory (or the Mooney theory with SEF U — (. ' (/ — 3)) leads for neo-Hookean materials to a true 
stress T = ^(A^ — j^-), or an engineering or nominal stress 


_T _E _l_ 
a ~ \ 1 ~ 3 (Al A 2 


where in terms of deformation u in the x direction we have, since deformations in the y and z directions are 
negligible in simple extension, 


A 


2 

1 


1 + 2c xx 

i+2 ^ 

ox 




or Ai — 1 + Note in this case the stress tensors T, S reduce to one nontrivial component T = Tn, 
<7 = Sn- Here E is a generalized modulus of elasticity and we note that, these formulations are restricted to 
Ai > 1 . 
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Figure 1: Rod with tip mass under tension. 


This can be used in the Timoshenko theory for longitudinal vibrations of a rubber bar with a tip mass 
to obtain 


pA c 


di 


da 


dt 2 
d 2 - 


■7 


dt dx 


= 0 0 < x < 




(1) 


F(t) + Mg, 


X—i 


where p is the mass density, F(t) is the applied external force, A c is the cross-sectional area, M is the tip 
mass, g is the gravitational constant, 7 is the air damping coefficient and 


a 


AcE <\ dr 5Al 

-(A 1 --)+A c C D — 


3 

a c e . 


-9 


du\ d 2 u 

d^ + c D ~dtdx 


( 2 ) 


is the internal (engineering) stress resultant, with g(^) = 1 + ^—(1+^) -2 . Here we have included a Kelvin- 
Voigt damping term \Cd is the Kelvin-Voigt damping coefficient) in the stress er as a first attempt to model 
damping. This leads to the nonlinear partial differential equation 


d 2 u du 

pAc W + 1 ~dt 


d_ 

dx 



+ A c Cn 



= 0 


0 < x < £ 


(3) 


d^u 


^t( d ^]+moAF\ 


dx 


dtdx J 


F(t) + Mg 


X—i 


for dynamic longitudinal displacements of a neo-Hookean material rod in extension. Since a series expansion 
°f 9 yields <j(£) = 3^ — 3 ^ 2 + 4 ^ 3 — • • •, this is readily seen, in the case of small displacements, to reduce to 
the usual longitudinal deformation equation for Hookean materials. 

In general, one does not expect the initial boundary value problem associated with (3) to have a classical 
(smooth) solution. Equations such as (3) are infinite dimensional, and must be discretized before the 
solution can be approximated. We have chosen to use a Galerkin method with linear splines for the spatial 
discretization. The second order equation (3) is then written as a first order system in time, and a stiff 
equation solver is used for the time integration. 

As reported in [3], the form of g for (3) obtained by either the neo-Hookean or simple Mooney SEF is 
inadequate in capturing the behavior of most filled elastomers. An essential task then is to determine a more 
general nonlinearity g with the aid of experiments and inverse problem techniques. 
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2.2 Approximation of nonlinear constitutive laws 

The neo-Hookean model (3) provides a natural example of a nonlinear PDE in the modeling of elastomers, 
but has only limited practical application since it is inadequate in describing most filled elastomers. Typically, 
one would employ equations such as (2) with a more general nonlinearity g which should be estimated from 
experiments. One does not expect such a general g to admit a SEF as a function of either the strain invariants 
or the extension ratios. Comparisons with SEF methods can be made by using the (approximate) SEF to 
derive the expected stress-strain relationship, and comparing results in the stress-strain plane. With this 
goal in mind, we now proceed to discuss the numerical estimation of g. 

2.2.1 Quasi-static experiments and inverse problems 

Although our ultimate goal is to use the results of dynamic experiments to determine the nonlinearity g, a 
reasonable first step is to estimate g using data from quasi-static pull tests. While this type of experiment 
cannot be used to study damping, it can be used to study the nonlinearity g, and more generally, the 
constitutive stress-strain law. 

With this in mind, we designed and implemented quasi-static pull tests on an Instron machine at Lord 
Corporation. A cylindrical rubber sample is secured vertically into the test machine, so that the lower end 
of the rod (x = 0) is fixed and the upper end of the rod (a: = £) is attached to a load cell and a horizontal 
crossbar (see Figure 2). A displacement pattern A (t) for the upper end of the rod is programmed into the 
machine, and the resulting loading force f(t) is recorded. In these experiments we initially used a constant 
rate of displacement at 5 inches per minute. 

Using the quasi-static load and displacement data at x = £, an inverse problem can be formulated to 
estimate the nonlinearity g. In general, this type of inverse problem is infinite dimensional both in state and 
parameter space, and hence finite dimensional approximations must be made. Preliminary results suggested 
a nonlinear but piecewise linear form (e.g., linear splines) might perform well in approximating the function 
g. This type of approximation for g was able to capture the load and displacement curves for unfilled natural 
rubber, re-confirming the need for a nonlinear g (see [3] for details). 

After obtaining satisfactory results for the quasi-static case, we may use the resulting form for g as an 
initial estimate in algorithms for determining g with dynamic data. We next describe these results. 

2.2.2 Dynamic experiments and inverse problems 

With the aid of scientists at Lord Corporation, we designed and carried out dynamic free-release and impulse 
hammer experiments with various samples of elastomers. These two types of experiments were chosen not 
only for the relative simplicity of the model for simple extension; it is well known that even in simple 
engineering materials (e.g., a spring-mass-dashpot system), it is easier to study dissipation of energy from 
free release or impulse response data than from other types of experiments such as cyclical tests (the almost 
universal dynamic tests currently reported in rubber research literature). 
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Linear fit to load cell cdata for unfilled rubber sample, small weiaht Linear fit to load cell cdata for unfilled rubber sample, small weiaht 



Figure 3: (left) Time domain approximation with a linear g, and (right) the FFT of the solution and the 
data. 


Four term fit to load cell data for unfilled rubber sample, small weight 
3|-,-,-r-,-.-.-.— 



Time (seconds) 



Frequency (hz) 


Figure 4: (left) Time domain approximation with a four-term piecewise linear g, and (right) the FFT of the 
solution and the data. 


In these experiments, a cylindrical elastomeric rod is suspended vertically so that the top end (a: = 0) is 
fixed, and a tip mass is attached at the lower end (a: = £) as in Figure 1. The frame at the bottom of the 
structure was used as an additional mass and as a housing for an accelerometer. An additional accelerometer 
is placed above the top of the elastomer to verify the clamped boundary condition at x = 0. A load cell is 
attached between the elastomer and the top accelerometer to provide force measurements at x = 0. 

For the free release experiment, the frame and rod were supported so that the rod itself was at its natural 
length (i.e., no compression or extension). The support was then removed, allowing the mass to fall freely. 
In the impulse response experiment, the rod and mass hung freely until equilibrium was reached, and then 
a hammer was used to excite the system. 

In both experiments described above, we used an unfilled natural rubber, natural rubber very lightly filled 
with carbon black, and highly filled natural rubber samples. In this section we present results pertaining 
to the unfilled sample. Using insight gained from the quasi-static inverse problem, we used the following 
inverse problem formulation to obtain the nonlinearity g: 

1 M du 

Minimize J(g) = - ^ |z; - A c a( , 0; g)) | 2 
i =1 

over g in some admissible class of functions Q. The observations z;, obtained from the experiments, are the 
force measurements at the fixed end, and u is the solution of (3) corresponding to g. 

Two forms for g were chosen for use in the inverse problem. First g was chosen as the best possible linear 
function, and this was followed by choosing the best g from a class Q of four-term piecewise linear splines. 
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Figures 3 and 4 demonstrate that the best linear g does not adequately capture the dynamics of the rubber, 
while a nonlinear function can give satisfactory results. 

These same experiments and inverse problems were repeated for medium and highly-filled rubber, but 
with less encouraging results. As the amount of filler is increased, the hysteretic behavior of the elastomer 
is more prevalent and must be included in the model. Since the highly-filled elastomers are more frequently 
used ill engineering applications than unfilled rubbers because of their increased strength and stiffness, an 
accurate model for the dynamic behavior of highly-filled elastomers is most desirable. This led to the need 
for developing a nonlinear, hysteretic constitutive law for filled elastomers. 

3 Nonlinear and hysteretic models, experiments and results 

The modeling of hysteresis in viscoelasticity is an often-studied but largely unresolved issue. Among the 
literature on modeling viscoelastic materials (e.g., [8], [10], [13], [24], [28]), two major types of stress-strain 
laws exist. One model is phenomenological, based on the mechanical behavior of the materials, and the other 
model stems from the microscopic behavior of fillers in the material. One common characteristic of existing 
constitutive models is that their formulations arose largely from experimental observations. 

The most fundamental of the phenomenological models is the Boltzmann integral model, which is formu¬ 
lated to capture dependence of the stress on the strain and/or strain rate. Experimental studies of elastomers 
have demonstrated a dependence of the stress on the strain and strain rate histories, suggesting that the 
Boltzmann integral formulation may well be appropriate here. One form of the Boltzmann integral model is 
given by 

a(t) = f Y(t — s) ^ ds, 

J— co ds 

where a(t) and s(t) are the stress and (infinitesimal) strain at time t respectively, and Y is known as the 
relaxation modulus function (see [7], [28] and the references therein). It can be shown that this type of 
constitutive relation is equivalent to certain internal variable models, such as those found in [15],[18], and 

[19]- .. 

Many of the existing models found in the literature have been verified for specific materials with individual 
experiments and computations. These results are difficult to generalize to other types of experiments and 
materials, however, since rubber materials are heavily dependent on various physical parameters. Our 
attempts to capture the dynamic behavior of elastomers include the use of a Boltzmann integral model, 
which we first developed in the quasi-static case before including it in the full dynamic model. 

3.1 Quasi-static hysteresis loops 

As seen in the previous section, a reasonable first step in developing a hysteretic constitutive law for (1) 
is to return the quasi-static case. Here we develop a model based on a Boltzmann integral term, with the 
inclusion of two types of nonlinearities. The form for our model was largely determined by experimental 
observations of quasi-static stress-strain curves for various samples of elastomers, and our inability to obtain 
useful fits to data with linear hysteresis models. 

One well-known characteristic of rubber-like materials is that they exhibit nested stress-strain curves, 
known as hysteresis loops (see Figure 5). This type of data was collected on the Instron machine (as in 
Section 2.2.1), and Includes a sequence of loading and unloading the sample with progressively smaller 
maximum strain levels (i.e., the outer-most loop is created first, then the inner loops follow in sequence). 

Our main approach to modeling is a pseudo-phenomenological one. We assume that the elastomer’s 
quasi-static behavior is prescribed by both an elastic response and a viscoelastic response. We choose a law 
of the following form: 

<r{t) = 9e(s{t)) + j Y(t - s)^-g v (£(s),£(s))ds, (4) 
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Figure 5: Nested hysteresis loops for highly-filled carbon black. 


where g e and g v are the elastic and viscoelastic response functions respectively, and Y is known as the 
memory kernel. Experiments and calculations in the quasi-static case led to the conclusion that elastomers 
have a finite memory, so that the stress depends on the strain and strain rate only for a history of length r. 
That is, the memory kernel Y obeys Y(t) 0 for all t > r, which implies that we may approximate (4) by 

c(0 = 9e{e{t)) + [ Y(t — s) -j-g v (£(«), i(s))ds. (5) 

Jt-r d‘ s 


The quasi-static stress-strain behavior in Figure 5 suggests that the elastomer possesses a different vis¬ 
coelastic response for loading (e > 0) than for unloading (e < 0). As a result, we choose a piecewise 
continuous form for the function g v : 


g v (£(s),£(s)) 


for e(«) > 0 

gvd{s{s)) for e(s) < 0 . 


We define the points t*, k — 1,. K to be the “turning points,” or the points in time for which i = 0. In 
our definition of g v , we do not require g v to be continuous at the turning points, so the derivative in (5) 
must be interpreted in the distributional sense. 

If we assume that the rubber is unloading for fj, < t < t^+i; where k is odd, then we can integrate by 
parts in (5): 


a(t) = 0) ‘ / Y(t - s)g v (s(s),i(s))ds+ Y(0)g vd (s(t)) 

Jt-r 

K 

+ "y]Y(t- t k ){-l) k+ 1 [g v i{e{tk)) - g V d{e(tk ))] 

k — 1 

where tx <t < t k+i, K odd. A similar expression holds for loading (t i< < t < tx+ i, K even). 

Different forms for g e , g v and Y were tested in comparison with experimental data, including the special 
cases of a linear g e and g v , with g v i = g v d- Our current model utilizes an exponential memory kernel of the 
form 

Y(t) = C,< 


and we have 




i = 0 


9v 


i x ) = E 


a*x 


i=0 


9vd(%') — ^ ^ b{X , 
i=0 




AP225 



Figure 6: Nested hysteresis loops for highly-filled carbon black: data and model predicted curves. 


where Ci, C% > 0 and K-i , (/,-, b . i = 0,1, 2 , 3 are real constants. 

The constants Ci,C 2 and Ek,(ik,bk are material-dependent parameters that must be determined us¬ 
ing parameter estimation techniques. This is implemented by formulating an inverse problem based on 
quasi-static experimental data. Using the Instron as before, we obtain displacement and force observations 
(A(ij), i = 1,..., N at the point x = l. These observations can be converted into infinitesimal strain 

i(ti) and stress <t(U) according to e — A/I and a = f/A c , where A c is the original cross-sectional area of 
the sample. 

We estimate the parameters Ci,C 2 , Ek,ak,bk, k = 0,1,2, 3 (which are collectively denoted as q ) by 
fitting this stress-strain data in the following inverse problem: 

1 N 

Minimize J (q) = - ^ 1^0 “ a ( t i 5 1) f ( 6 ) 

i =1 

over q in some admissible parameter set Q, where tr(ti',q ) denotes the quantity obtained by inserting the 
data s(ti) into (5) with the parameters q. 

We carried out the computations in MATLAB and FORTRAN, using both BFGS-type and Nelder- 
Mead optimization routines. This inverse problem was solved using data from several types of elastomers, 
including silica-filled silicone and lightly, medium and highly-filled carbon black. Here we present results for 
the highly-filled carbon black sample (the type of elastomer most important to industry). 

Our main goal in solving an inverse problem is to find parameters q* that will best predict a set of nested 
hysteresis loops like those found in Figure 5. If possible, we would like to use data from at most one or 
two of the loops in the inverse problem itself, yet still obtain an accurate prediction of the entire set. After 
experimenting with different data sets, we chose the data from the 100% and 90% strain loops (the two 
outer-most loops in Figure 5) for use in the inverse problem. That is, we let {e(U), cr(f,;)}, i = 1,.. ., TV in ( 6 ) 
be the strain and stress data from the 100% and 90% strain loops, and we solved the inverse problem for ( 6 ) 
to obtain an optimal parameter set q*. We then used these parameters in (5) to simulate stress-strain curves 
and compare with data for the inner loops. The resulting curves are shown in Figure 6 , with a relative error 
of 2 . 2 % between model prediction and data. 

Equivalent results for the other elastomer samples were similarly encouraging, leading us to the conclusion 
that the constitutive relation (5) provides satisfactory predictions of quasi-static stress-strain data for several 
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types of elastomers. Since our ultimate goal is to accurately predict the dynamic behavior of elastomers, we 
return to the full dynamic case and attempt to incorporate (5) into the dynamic model (1). 

3.2 A dynamic model with hysteresis 

Up to this point, none of the constitutive models used in the dynamic case were able to capture the hysteretic 
behavior of elastomers. After developing the stress-strain law (5) with quasi-static data which successfully 
predicts hysteresis in several types of elastomers, we are ready to proceed in including (5) in the full dynamic 
model. This however leads to some immediate theoretical questions. 

3.2.1 Theoretical issues 

By combining the constitutive relation (5) with (1), we obtain the following dynamic model for tx < t < 
tx+i, K odd (this is for unloading — a similar model holds for loading): 





du 

dx 


m(0, x ) 
u t (0, x) 
u(t, 0 ) 
u(t, X ) 


d_ 

dx 


A, 


du 

dx 


. _ ,du d 2 u 

' cg ^^ + cCd dtdx 


■ A C Y ( 0 )g v d{-g^) + A c J 


4 ■ du 

Y(t - s)g v ( — )ds 


K 




*+l 


k — 1 


, du ,du, 

9vi ~dx ~ 9vd cfx ^ ■ 


, 0 < x < 


9e{ W {t)) + Cd U- x + nO) 9vd (^(t)) + J y(t - s)g v ( d £(s))ds 


(7) 

( 8 ) 


K 

k -1 


cd u /d u s w 

9 v i\~g^V'k)) — 9vd\~Q^\tk)) 


X—i 


m 


(9) 

( 10 ) 


f 0 


( 11 ) 


0 

0 

<pi, t < 0 . 


( 12 ) 

(13) 

(14) 


Here we assume that the elastomer, with the usual fixed end at x = 0, is at rest at time t = 0 with 
deformation ipo and memory <p\. In addition, we include an (internal) Kelvin-Voigt damping term . h't iiHsti 
( Cd > 0) in the model. 

At this stage, it is unclear whether the above PDE is in fact well-posed. That is, can we even guarantee 
the existence of a solution for (7) - (14)? This issue is fully addressed in [4]; we present an outline of the 
results here. 

Consider the Hilbert spaces V = 0,1) = {<j> : <j> £ i7 1 (0, l), <j>(0) = 0} and H = i 2 (0,7), so that we 

have the Gelfand triple V H V*. We denote (•, •) as the inner product in H , while {•, f)v*,V stands 
for the usual duality product. Let £[ a ,h\ denote the space of functions u : [a, b\ —>■ H such that. 

£[a,b] = {w : U e L 2 (a, b; V), u t 6 L 2 (a, b; V)}. 

Motivated by the model (7) - (14) we consider the system 

d fdu du .du. ■ .du. , 

««■ - c °- y x (^ + + /_ r y <* - 

+ E y ('-**)(-l)‘ +1 b.i(|j(<*))-9.4$)(!*))]) =m in V* (15) 

k — 1 / 
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S2 

ot>. 

O 

II 

o 

(16) 

£ 

II 

-6 

o 

(17) 

«j(a, ai) = 0 

(18) 

£ 

II 

C-+. 

A 

(19) 


for a <tK <t < tx+i = b. 

We say that u G £[ a j] is a weak solution of (15)-(19) with po £ h and p\ G L 2 (a — r , a,; H) if it satisfies 
/ [ - { u t ,<p t ) + {u x ,p x } + C D {u xt , + {g e (u x ),p x ) + {Y ((%, (v x . ii x ) : <f x ) 

J a 

K 

+(£h{t - tk)Y (t - 4)(-l) i+1 [ff„i( M ®(^)) “ fiw(M^))]>^) 

& = 1 

+{ / Y{t- s)g v {u x , u x )ds, <p x )]dr + p) = / {F,p}v,vdT (20) 

Jt—v Ja 


for any a < < t < tjf+i = b and ^ G £[ a ,b] as well as the initial condition 

u(a, x) = tpo. 


( 21 ) 


Here we assume that g e , g v i and g v d are continuous nonlinear mappings of real gradient type, so that 
there exist continuous Frechet differentiable nonlinear functionals which obey a certain type of boundedness. 
Moreover, the nonlinear functions g e , g v i and g v( j are bounded in a particular sense, and their Frechet 
derivatives are bounded in C(H, H). We also assume that g e , g v i and g v d obey a monotonicity condition, and 
that Y is a smooth and bounded function. For details and more complete statements of these assumptions, 
see [4]. 

It is shown in [4] that under the precise assumptions detailed there, there exists a global weak solution of 
(15)-(19) for any ipo G V, G L 2 (a — r, a; V). We first proved that under these assumptions a weak solution 
G £[o,«i] exists for any (po G V and ip i G L 2 (—r , 0; V ), i.e., we used a = to = 0, b = ti in our definition 
of the weak solution. We showed that t/ 1 ' is smooth enough so that u {1 \ti), (<■) exist. Next we 

suppose that a unique weak solution u( K \t) exists on [0,1 a'] and we consider the interval [IajIa'+i]- Then 
we have a similar system as above except that now we pick up a jump term and po and p\ are modified, 
i.e., ipo = u( K \tx) and 


j pi if t < 0 
\ if 0 < t < t K . 


We again show that a unique weak solution u exists on this interval with the necessary smoothness and 



if t <tx 

if t K < t < t K +i- 


is a weak solution on [ 0 ,tA+i]- 

To show the existence of the weak solution uW on [0,fi], we first give a priori estimates, then introduce 
Galerkin approximations to the weak solution and justify taking the limit. We used the Minty-Browder 
technique in showing that the limit of the Galerkin approximates is a weak solution. A complete statement 
of the theorem and its proof can be found in [4]. 

3.2.2 Numerical results 

After addressing the theoretical issues associated with the dynamic, hysteretic model, we tested the effec¬ 
tiveness of this model on our previously obtained dynamic elastomer data. We used the free release dynamic 
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Figure 7: Dynamic results for lightly-filled carbon black. 


experiments with a 3 lb. tip mass attached to the free end of a lightly filled carbon black sample. To 
model this particular sample we use a cubic polynomial for </ e (e) = die + a 2 e 2 + (as in the quasi-static 
case), and we suppose that g v is linear and does not depend on £, i.e., g v — g v i = g vl j- As before, we used 
V'(/) = C 2 e~ Clt , where C\,C 2 > 0. 

We note that although our dynamic model (7) contains Kelvin-Voigt damping, we do not include it in our 
computational model. This damping term is important for the theoretical result, but in our experiments and 
computations we have determined that the hysteresis portion of the model can provide adequate damping 
to capture the dynamic behavior. We therefore can take the Kelvin-Voigt coefficient as Cp 0, and thus 
reconcile the theoretical and numerical results. 

Using the load cell data from the experiment, we set up a parameter identification problem to estimate 
the parameters p,ai,a 2 ,as,Ci,C 2 (which we will denote collectively as q). That is, we wish to minimize 

a 2 

.Oil, .. 

-A-cO "(^^7 (A 5 (?)) J 

where z;, i = 1,.. ,,M are observations from the load cell. Moreover, u is a solution to (7) with parameters 
9- 

In our computations, we used linear splines for spatial discretization, while piecewise constant elements 
were used in the hysteretic discretization. (More details on the computational technique for dealing with 
the integral term can be found in [1],[2].) We used MATLAB optimization routines for the inverse problem. 
The computed result shows very good agreement with the collected data, with a relative error of 3% (see 
Figure 7). Current efforts involve use of these models and techniques with experimental data from highly 
filled samples. 

4 Conclusion 

In the above presentation we have outlined the collaborative effort between applied mathematicians in the 
CRSC and scientists at Lord Corporation, focused on the modeling of elastomers. We have made significant 
progress to date towards the final goal of developing a high-quality model for the dynamic behavior of 
elastomers. 


j(t) = J2 
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The initial phase of the project was an exploration of nonlinear constitutive laws. We considered several 
types of nonlinear stress-strain laws, including many found in the existing literature. Working with both 
the quasi-static and the dynamic cases, our experiments and computations suggested that neither the neo- 
Hookean nor a linear constitutive relation are adequate for capturing the dynamic behavior of elastomers, 
even for unfilled natural rubber. For very lightly-filled and unfilled rubber, it appears that a nonlinear 
stress-strain law is sufficient, while for more highly-filled rubbers it is clear that hysteresis must be included 
in the constitutive model. 

The next phase of our collaboration involved incorporating hysteresis into our quasi-static and dynamic 
models. We refined our stress-strain law first using quasi-static data, until we obtained a model that could 
accurately predict the quasi-static behavior of several types of elastomers. We then returned to the dynamic 
case to incorporate the hysteretic constitutive relation in our dynamic model. After formulating the full 
dynamic PDE model with hysteresis, it became clear that the question of the well-posedness of the model 
needed to be resolved in a formal manner. This led to the development of the weak form of the PDE and a 
proof of the existence of solutions. Once these theoretical questions were addressed, we tested the accuracy 
of this model in the dynamic case by using inverse problem techniques. As shown in the previous section, 
this model produces satisfactory results for lightly-filled rubber. Current work includes the testing of our 
dynamic model on more highly-filled rubbers, which exhibit more significant hysteresis. 

Future work in this collaboration includes a study of the mechanical behavior of elastomers in shear. 
Since the overall goal of the project is to accurately predict the dynamic behavior of elastomers, we need 
to understand both extension and shear. As with extension, we will work with both the quasi-static and 
dynamic paradigms to develop a model for shear, and we will test the accuracy of our model by gathering 
shear data for various types of elastomers. 
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